import pandas as pd

# Load your Excel file
file_path = r"rainfall.xlsx"  # Replace with your actual file path
data = pd.read_excel(file_path)

# Save the result to a new Excel file
output_file = "average_monthly_totals.xlsx"

# Convert the 'Date' column to datetime
data['Date'] = pd.to_datetime(data['Date'])

# Extract the Month and Year from the 'Date' column
data['Month'] = data['Date'].dt.month
data['Year'] = data['Date'].dt.year

# Reshape the data from wide to long format
melted_data = data.melt(id_vars=['Date', 'Month', 'Year'], 
                        var_name='Station_ID', 
                        value_name='Rainfall')

# Convert 'Rainfall' to numeric, coercing errors to NaN
melted_data['Rainfall'] = pd.to_numeric(melted_data['Rainfall'], errors='coerce')

# Drop rows with missing rainfall data
melted_data = melted_data.dropna(subset=['Rainfall'])

# Step 1: Calculate total monthly rainfall for each station for each year
monthly_totals = (
    melted_data.groupby(['Station_ID', 'Year', 'Month'])
    .agg(Total_Monthly_Rainfall=('Rainfall', 'sum'))
    .reset_index()
)

# Step 2: Calculate the average of total monthly rainfall across all years for each station
average_monthly_totals = (
    monthly_totals.groupby(['Station_ID', 'Month'])
    .agg(Average_Monthly_Rainfall=('Total_Monthly_Rainfall', 'mean'))
    .reset_index()
)

#Saving file
average_monthly_totals.to_excel(output_file, index=False)

print(f"Average monthly rainfall totals saved to {output_file}.")

import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.lines import Line2D

# Load Excel file
file_path = r"average_monthly.xlsx"

# Load the shapefile (replace the file path with your shapefile path)
shp_file = r"gadm41_IDN_1.shp" 

# Function to convert decimal degrees to DMS format with N/S/E/W suffix
def decimal_to_dms(decimal_degree, is_latitude=False):
    # Get the degree part
    degrees = int(decimal_degree)
    # Get the minutes part
    minutes = int((decimal_degree - degrees) * 60)
    # Get the seconds part
    seconds = (decimal_degree - degrees - minutes / 60) * 3600
    # Round to 2 decimal places for seconds
    seconds = round(seconds, 0)
    
    # Append N/S for latitude or E/W for longitude
    if is_latitude:
        direction = 'N' if degrees >= 0 else 'S'
    else:
        direction = 'E' if degrees >= 0 else 'W'
    
    return f"{abs(degrees)}°{abs(minutes)}'{abs(seconds)}\" {direction}"

# Data load from Excel
avg_data = pd.read_excel(file_path, sheet_name='Avg') # Make sure the sheet name is correct
coord_data = pd.read_excel(file_path, sheet_name='Coord') # Make sure the sheet name is correct

# Merge the data
merged_data = pd.merge(avg_data, coord_data, on='ID') # Make sure the column name is correct

# Create geometry column
merged_data['geometry'] = merged_data.apply(lambda row: Point(row['lon'], row['lat']), axis=1)
gdf = gpd.GeoDataFrame(merged_data, geometry='geometry')

# Get min and max rainfall across all stations 
min_rainfall = gdf['Avg Rainfall'].min() # Make sure the column name is correct
max_rainfall = gdf['Avg Rainfall'].max() # Make sure the column name is correct

# Define custom intervals based on your specifications
intervals = [min_rainfall, 100, 150, 200, 250, 300, 350, 400, max_rainfall] # Customize as needed

# Function to assign color based on rainfall amount
def get_color(rainfall):
    if rainfall <= intervals[1]:
        return 'yellow'
    elif intervals[1] < rainfall <= intervals[2]:
        return 'red'
    elif intervals[2] < rainfall <= intervals[3]:
        return 'green'
    elif intervals[3] < rainfall <= intervals[4]:
        return 'pink'
    elif intervals[4] < rainfall <= intervals[5]:
        return 'blue'
    elif intervals[5] < rainfall <= intervals[6]:
        return 'purple'
    elif intervals[6] < rainfall <= intervals[7]:
        return 'orange'
    else:
        return 'brown'

# Set SHP as basemap
base_map = gpd.read_file(shp_file)

# Map the marker sizes based on rainfall intervals
def get_marker_size(rainfall):
    """
    Adjust the size of the markers based on the rainfall intervals.
    The size is scaled according to the defined intervals.
    """
    if rainfall <= intervals[1]:
        return 4  # Small size for the lowest range
    elif intervals[1] < rainfall <= intervals[2]:
        return 8
    elif intervals[2] < rainfall <= intervals[3]:
        return 12
    elif intervals[3] < rainfall <= intervals[4]:
        return 16
    elif intervals[4] < rainfall <= intervals[5]:
        return 24
    elif intervals[5] < rainfall <= intervals[6]:
        return 32
    elif intervals[6] < rainfall <= intervals[7]:
        return 40
    else:
        return 46  # Largest size for the highest rainfall range

#'Avg Rainfall' column which contains the rainfall data
fig, ax = plt.subplots(figsize=(12, 12))

# Plot the high-resolution base map from SHP file
base_map.plot(ax=ax, color='white', edgecolor='black', linewidth=0.5)

# Apply marker size based on rainfall and plot the data points
gdf.plot(ax=ax, markersize=gdf['Avg Rainfall'].apply(get_marker_size),
         marker='o', color=gdf['Avg Rainfall'].apply(get_color), 
         edgecolor='black', linewidth=0.5,  
         label='Rainfall')

# Expand the focus to a larger region around Indonesia, customize as needed
ax.set_xlim(93, 143)   
ax.set_ylim(-13, 9)    

# Create the legend and place it outside the plot area
legend_labels = [f'{round(intervals[0], 2)} - {round(intervals[1], 2)}', 
                 f'{round(intervals[1], 2)} - {round(intervals[2], 2)}',
                 f'{round(intervals[2], 2)} - {round(intervals[3], 2)}',
                 f'{round(intervals[3], 2)} - {round(intervals[4], 2)}',
                 f'{round(intervals[4], 2)} - {round(intervals[5], 2)}',
                 f'{round(intervals[5], 2)} - {round(intervals[6], 2)}',
                 f'{round(intervals[6], 2)} - {round(intervals[7], 2)}',
                 f'{round(intervals[7], 2)} - {round(intervals[8], 2)}']  

legend_colors = ['yellow', 'red', 'green', 'pink', 'blue', 'purple', 'orange', 'brown']

# Create custom legend handles with marker size scaling for legend
legend_marker_sizes = [4, 5, 6, 7, 8, 9, 9.5, 11] # Customize as needed

# Adjust the maximum size of the legend markers to avoid exceeding the legend box
max_legend_marker_size = 11 
handles = [Line2D([0], [0], marker='o', color='w', markerfacecolor=color, markersize=min(size, max_legend_marker_size), 
                  markeredgecolor='black', markeredgewidth=0.5) 
           for color, size in zip(legend_colors, legend_marker_sizes)]

# Position the legend outside the frame using bbox_to_anchor
plt.legend(handles=handles, labels=legend_labels, loc='upper left', bbox_to_anchor=(1, 1), title="BMKG (mm/month)", fontsize=9) # Adjust title name based on data

# Add title and show plot, customize as needed
plt.title("Average Monthly Rainfall Indonesia")

# Set the coordinate labels to DMS format with N/S/E/W
xticks = ax.get_xticks()
yticks = ax.get_yticks()

ax.set_xticklabels([decimal_to_dms(x, is_latitude=False) for x in xticks])
ax.set_yticklabels([decimal_to_dms(y, is_latitude=True) for y in yticks])

# Adjust layout to fit the plot and legend
plt.tight_layout()
plt.show()

import pandas as pd
import geopandas as gpd
import numpy as np

from shapely.geometry import Point
from shapely.geometry import mapping
from scipy.spatial import cKDTree
from scipy.interpolate import griddata

import matplotlib.pyplot as plt
from matplotlib.colors import Normalize

import rasterio
import rasterio.mask
from rasterio.features import geometry_mask
from rasterio.transform import from_origin
from rasterio.mask import mask

# Load Excel data
file_path = r"average_monthly.xlsx"
avg_data = pd.read_excel(file_path, sheet_name='Avg')
coord_data = pd.read_excel(file_path, sheet_name='Coord')

# Load shapefile for the base map
shp_file = r"gadm41_IDN_1.shp"

# Merge the data
merged_data = pd.merge(avg_data, coord_data, on='ID')   # Make sure the column name is correct

# Extract coordinates and values
coordinates = merged_data[['lon', 'lat']].values        # Make sure the column name is correct
values = merged_data['Avg Rainfall'].values             # Make sure the column name is correct

# Define grid for interpolation
x_min, x_max = coordinates[:, 0].min() - 1, coordinates[:, 0].max() + 1
y_min, y_max = coordinates[:, 1].min() - 1, coordinates[:, 1].max() + 1
grid_x, grid_y = np.meshgrid(
    np.linspace(x_min, x_max, 500),  # Adjust resolution here
    np.linspace(y_min, y_max, 500),  # Adjust resolution here
)

# Define the IDW function
def idw(x, y, values, xi, yi, power=2):
    """
    Perform IDW interpolation.
    
    :param x, y: Coordinates of known points.
    :param values: Known values at the points.
    :param xi, yi: Grid coordinates for interpolation.
    :param power: Weighting power.
    :return: Interpolated values on the grid.
    """
    tree = cKDTree(np.c_[x, y]) 
    distances, indices = tree.query(np.c_[xi.ravel(), yi.ravel()], k=len(x))
    weights = 1 / (distances ** power)
    weights[np.isinf(weights)] = 1e12  # Handle division by zero for very small distances
    weights /= weights.sum(axis=1, keepdims=True)  # Normalize weights along rows
    zi = np.sum(weights * values[indices], axis=1)
    return zi.reshape(xi.shape)

# Perform IDW interpolation
interpolated_values = idw(coordinates[:, 0], coordinates[:, 1], values, grid_x, grid_y)

# Set SHP as basemap
base_map = gpd.read_file(shp_file)

# Function to convert decimal degrees to DMS format
def decimal_to_dms(degrees):
    # Determine the direction (N/S or E/W)
    direction = 'N' if degrees >= 0 else 'S' if degrees < 0 else 'E' if degrees >= 0 else 'W'
    
    degrees = abs(degrees)  # Work with absolute value
    d = int(degrees)  # Degrees
    m = int((degrees - d) * 60)  # Minutes
    s = int((degrees - d - m / 60) * 3600)  # Seconds (converted to integer)
    
    return f"{d}°{m}'{s}\"{direction}"

# Set map boundaries and labels with DMS format
def format_ticks(ax):
    # Convert the ticks to DMS format
    x_ticks = ax.get_xticks()
    y_ticks = ax.get_yticks()
    
    # Convert to DMS for longitude (x) and latitude (y)
    x_labels = [decimal_to_dms(tick) for tick in x_ticks]
    y_labels = [decimal_to_dms(tick) for tick in y_ticks]
    
    # Set the new labels
    ax.set_xticklabels(x_labels, fontsize=7)
    ax.set_yticklabels(y_labels, fontsize=7)

# Plot the IDW raster
fig, ax = plt.subplots(figsize=(12, 12))
contour = ax.contourf(grid_x, grid_y, interpolated_values, levels=20, cmap='RdYlGn_r', alpha=1)

# Add the base map
base_map.plot(ax=ax, color='none', edgecolor='black', linewidth=0.5)

# Add station points for reference
plt.scatter(coordinates[:, 0], coordinates[:, 1], c=values, cmap='Reds', edgecolor='black', s=10)

# Add colorbar with smaller size
cbar = plt.colorbar(contour, ax=ax, fraction=0.019, pad=0.04)  # Adjust fraction and pad
cbar.set_label('IDW BMKG (mm/month)', fontsize=9)

# Set map boundaries and labels with DMS formatting
ax.set_xlim(x_min, x_max)
ax.set_ylim(y_min, y_max)
format_ticks(ax) 

plt.title("IDW Interpolated for Average Precipitation")
plt.tight_layout()
plt.show()